Progetto

Premesse

Il movente

Lo studio nasce dalla curiosità di una studentessa di toccare con mano alcune delle misure macroeconomiche italiane più importanti, spesso menzionate da economisti e media.

I dati

Tutti i dati analizzati nel seguente studio sono stati scaricati dal portale FRED, “Federal Reserve Economic Data”, un vasto archivio di dati gestito dalla Federal Reserve Bank of St. Louis.
I dati scelti sono tutti rappresentati da serie storiche, tutte a frequenza trimestrale e non stagionalmente adattate.
La finestra temporale prescelta è quella che va dal primo trimestre del 1996 al quarto trimestre del 2022 e i dati sono tutti riferiti al contesto economico italiano.

Gli strumenti

Gli strumenti statistici utilizzati sono quelli appresi durante il corso di Laboratorio di Statistica, incentrato sull’analisi delle serie storiche e sull’interpretazione dei relativi risultati.

I packages R necessari per l’analisi che si intende condurre sono i seguenti.

library(forecast)
library(quantmod)
library(GGally)
library(tidyverse)
library(dygraphs)

Real Gross Domestic Product

Il Gross Domestic Product (GDP), che in italiano traduciamo con Prodotto Interno Lordo (PIL), è una misura del valore di tutti i beni e servizi finali prodotti nel Paese in un determinato periodo di tempo. La misura Real GDP tiene conto dell’effetto dell’inflazione.
Il Real GDP riesce a fornire una misura accurata della crescita economica dell’Italia e può essere utilizzato per confrontare l’andamento dell’economia italiana con quello di altri Paesi o con il passato.

Nota: Ci si riferirà a questa serie storica sempre con l’acronimo GDP.

gdp <- ts(getSymbols("NGDPRNSAXDCITQ", src="FRED", auto.assign = FALSE), start=1996, frequency=4) 

Fonte: International Monetary Fund, Real Gross Domestic Product for Italy [NGDPRNSAXDCITQ], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/NGDPRNSAXDCITQ, May 11, 2023.

dyRangeSelector(dygraph(gdp, 
                        main="Real Gross Domestic Product",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

La prima informazione che è possibile ricavare dai dati e dalla loro visualizzazione riguarda l’andamento medio di lungo periodo della serie: il GDP Italiano ha avuto un andamento medio crescente fino alla metà del 2008, per poi mostrare un rapido calo fino agli inizi del 2009. Questa fase coincide con la crisi finanziaria del 2007-2008. Dopo una fase di ripresa, durata fino alla metà del 2011, il trend è tornato ad essere decresente fino agli inizi del 2014, per poi riprendere in salita fino alla fine del 2019. Nel 2020 l’emergenza Covid-19 ha avuto effetti drammatici sull’economia di ogni Paese e il GDP italiano ha toccato rapidamente, nella metà del 2020, i minimi storici. Tuttavia, la ripresa è stata rapidissima, mostrando un trend crescente ad alta pendenza.

Pulizia dei dati

Valori mancanti (missing values) o anomali (outliers) potrebbero avere effetti sui modelli che si andranno a costruire e sui relativi risultati. Risulta pertanto indispensabile verificare l’eventuale presenza di tali valori.

Missing values

anyNA(gdp)
## [1] FALSE

Outliers

tsoutliers(gdp)
## $index
## [1] 95 96 98
## 
## $replacements
## [1] 399055.6 397351.1 386198.6

La serie storica oggetto di studio non presenta missing values, pertanto non sarà necessario gestirli con metodi di imputazione; tuttavia, presenta outliers. I valori anomali presenti nella serie sono frutto degli effetti economici dell’emergenza Covid-19.

La funzione tsclean() utilizza la decomposizione robusta STL per la stima dei missing values e la sostituzione degli outliers. Possiamo utilizzarla per creare una serie storica in cui viene applicata la sostituzione degli outliers, in modo da poter confrontare i risultati dei metodi che verranno applicati sulle due serie.

gdp_clean <- tsclean(gdp)
autoplot(gdp_clean, series="Clean") +
  autolayer(gdp, series="Original") +
  ggtitle("Original series VS Clean series") +
  xlab("Year") +
  ylab("Currency (€)") +
  guides(colour=guide_legend(title="Series"))

Stagionalità

Si prosegue con la ricerca di un’eventuale stagionalità all’interno della serie storica del GDP.

ggsubseriesplot(gdp) +
  ylab("Currency (€)") + 
  ggtitle("Seasonal subseries plot: Real Gross Domestic Product")

ggseasonplot(gdp, year.labels = TRUE, year.labels.left = TRUE) +
  ylab("Currency (€)") + 
  xlab("Quarter")

  ggtitle("Seasonal plot: Real Gross Domestic")
## $title
## [1] "Seasonal plot: Real Gross Domestic"
## 
## attr(,"class")
## [1] "labels"

Come si poteva prevedere, il GDP può essere influenzato dalla stagionalità in quanto le attività economiche possono variare in modo regolare e prevedibile a seconda delle stagioni: il GDP mostra un andamento crescente il primo e il terzo trimestre, al contrario di ciò che accade il secondo e il quarto. Tuttavia, sembra che negli ultimi anni (orientativamente dal 2014) questo andamente sia stato meno netto e gli anni 2020 e 2021 hanno avuto andamenti anomali per le stesse considerazioni sopracitate legate all’emergenza Covid-19.

Decomposizione

É possibile decomporre la serie storica GDP isolando le diverse sue componenti: trend-ciclo, stagionalità e residuale. Si procede con la decomposizione STL (Seasonal and Trend decomposition using Loess), in quando è la più robusta ai valori anomali.

stl_fit <- stl(ts(as.vector(gdp), start=1996, frequency = 4), s.window = "periodic", robust=TRUE)
autoplot(stl_fit)

La decomposizione applicata ai dati originali appare fortemente influenzata dai valori anomali, nonostante la decomposizione STL sia la più robusta. Pertato, si applica la decomposizione alla serie gdp_clean.

stl_fit2 <- stl(ts(as.vector(gdp_clean), start=1996, frequency = 4), s.window = "periodic", robust=TRUE)
autoplot(stl_fit2)

La componente principale risulta essere quella di trend.

autoplot(gdp) +
  autolayer(stl_fit$time.series[,2], series="original data") +
  autolayer(stl_fit2$time.series[,2], series="clean data") +
  ylab("Currency (€)") +
  ggtitle("STL decomposition") +
  guides(colour=guide_legend(title="STL on"))

Possiamo, ora che sono state separate le componenti di trend e stagionalità, calcolarne la forza.

var_R <- var(stl_fit2$time.series[,3])
var_RT <- var(stl_fit2$time.series[,3] + stl_fit2$time.series[,2])
var_RS <- var(stl_fit2$time.series[,3] + stl_fit2$time.series[,1])

Forza del trend:

\[ F_t=max(0;1-\frac{Var(R_t)}{Var(R_t+T_t)}) \]

max(0,1-(var_R/var_RT))
## [1] 0.9661618

Forza della stagionalità:

\[ F_s=max(0;1-\frac{Var(R_t)}{Var(R_t+S_t)}) \]

max(0,1-(var_R/var_RS))
## [1] 0.8836213

Entrambi gli indici sono compresi tra 0 e 1. Essendo entrambi prossimi a 1, \(F_t=0.97\) e \(F_s=0.88\), entrambe le componenti sono forti, ma il trend risulta essere la componente principale.

Autocorrelazione

Per misurare la forza del legame lineare tra i valori della serie storica su due istanti di tempo diversi, t e t-h, dove h è detto ritardo o lag, si procede con lo studio dell’autocorrelzione della serie del GDP italiano.
L’autocorrelazione è data da

\[ \rho(h)= \frac{\gamma(h)}{\gamma(0)} \] dove

\[ \gamma(h) = \frac{1}{T} \sum_{t=1}^{T-h} (x_t - \bar{x})(x_{t+h} - \bar{x}) \]

e \(\gamma(0)\) si riduce a

\[ \gamma(0) = \frac{1}{T} \sum_{t=1}^{T} (x_t - \bar{x})^2 \] ovvero alla nozione di varianza.

L’insieme di tutti i valori di \(\rho(h)\) è detto funzione di autocorrelazione. Si tratta di una funzione simmetrica intorno allo zero.

ggAcf(gdp, plot=FALSE)
## 
## Autocorrelations of series 'gdp', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.551  0.697  0.393  0.658  0.240  0.403  0.137  0.428  0.082  0.283 
##     11     12     13     14     15     16     17     18     19     20 
##  0.061  0.365  0.030  0.213 -0.045  0.236 -0.075  0.099 -0.137  0.141
ggAcf(gdp) 

Dai valori delle autocorrelazioni, calcolati per diversi ritardi, risulta evidente che le autocorrelazioni di ritardo \(h=2\) sono quelle dal valore maggiore, coerentemente con le considerazioni fatte precedentemente circa la stagionalità della serie. Ciò si riflette anche nel grafico che mostra regolarità nel comportamento delle autocorrelazioni.
I valori delle autocorrelazioni permettono di rilevare la presenza di trend, in quanto inizialmente sono tutte positive e decacono a zero lentamente, e del relativo cambio di pendenza, in quanto c’è un cambio di direzione per i valori delle autocorrelazioni.
Nonostante i valori delle autocorrelazioni nella prima parte del grafico ricadano al di fuori dell’intervallo di confidenza (fissato al 95%), mostrando delle autocorrelazioni statisticamente significative, queste diventano statisticamente nulle nella parte finale, per valori elevati di h.

Previsioni

I dati di una serie storica, in questo caso del GDP italiano, possono essere utilizzati per fare previsione sul valore futuro che questa variabile assumerà. Chiaramente l’intervallo di confidenza diventa via via più ampio all’allontanarci dall’ultima osservazione a disposizione, rendendo le previsioni a lungo termine più incerte.

La presenza di valori anomali influisce notevolmente sulle capacità predittive dei modelli applicati. I modelli verranno addestrati sulle osservazioni dal 1996 al 2018 e testati sulle osservazioni dal 2019 al 2022 della serie storica trasformata gdp_clean per ridurre l’effetto dell’emergenza Covid sui dati.

Modelli benchmark

Di seguito verranno implementati i più semplici metodi di previsione, seppur con poche pretese conoscendo le caratteristiche della serie:

  • Average method

\[ y_{T+h|T} = \frac{\sum_{t=1}^{T} y_t}{T} \ \]

  • Naïve method

\[ y_{T+h|T} = y_T \]

  • Seasonal naïve method

\[ y_{T+h|T} = y_{T+h-m(k+1)} \]

  • Drift method

\[ y_{T+h|T} = y_T + h \frac{{y_T-y_1}}{T-1} \]

train <- ts(gdp_clean[1:92], start=1996, frequency=4)
test <- window(gdp_clean, start=2019, end=c(2022,4))
mean_fit <- meanf(train, 16)
naive_fit <- naive(train, 16)
snaive_fit <- snaive(train, h=16)
drift_fit <- rwf(train, 16, drift=TRUE)
autoplot(gdp_clean) + 
  autolayer(mean_fit, series="Mean", PI=FALSE) + 
  autolayer(naive_fit, series="Naïve", PI=FALSE) +
  autolayer(snaive_fit, series="Seasonal Naïve", PI=FALSE) + 
  autolayer(drift_fit, series="Drift", PI=FALSE) +
  xlab("Year") + 
  ylab("Currency (€)") +
  ggtitle("Real Gross Domestic Product") +
  guides(colour=guide_legend(title="Forecast"))

rbind(
"Average" = accuracy(mean_fit, test)[2,c(2,3,5,6)],
"Naïve" = accuracy(naive_fit, test)[2,c(2,3,5,6)],
"Seasonal Naïve" = accuracy(snaive_fit, test)[2,c(2,3,5,6)],
"Drift" = accuracy(drift_fit, test)[2,c(2,3,5,6)])
##                    RMSE      MAE     MAPE     MASE
## Average        13891.99 11374.26 2.863320 1.672807
## Naïve          21352.64 16987.33 4.393137 2.498318
## Seasonal Naïve 13754.56 10821.43 2.785163 1.591503
## Drift          25741.42 22480.21 5.771592 3.306153

Chiaramente nessuno di questi metodi di previsione risulta efficace, ma tra tutti il seasonal naïve method è l’unico che riesce a cogliere la stuttura stagionale della serie, e mostra il più basso tasso di errore. Il seasonal naïve method è incapace di cogliere la componente di trend presente nei dati, componente colta dal drift method. Il metodo che garantisce un’accuratezza previsiva prossima a quella garantita dal seasonal naïve method, è l’average method, pur essendo il modello più semplice.

No Covid Data

Nonostante i modelli siano stati testati sulla serie gdp_clean, nessun metodo di previsione addestrato sulla serie fino al 2019 sarebbe stato in grado di prevedere quando accaduto nel 2020 al GDP italiano. Può essere utile valutare l’accuratezza di questi metodi addestrandoli sulle osservazioni dal 1996 al 2016 e testandoli sulle osservazioni dal 2017 al 2019 della serie storica originale gdp.

gdp_nc <- ts(gdp[1:96], start=1996, frequency=4)
train_nc <- ts(gdp_nc[1:84], start=1996, frequency=4)
test_nc <- window(gdp_nc, start=2017)
mean_fit_nc <- meanf(train_nc, 12)
naive_fit_nc <- naive(train_nc, 12)
snaive_fit_nc <- snaive(train_nc, h=12)
drift_fit_nc <- rwf(train_nc, 12, drift=TRUE)
autoplot(gdp_nc) + 
  autolayer(mean_fit_nc, series="Mean", PI=FALSE) + 
  autolayer(naive_fit_nc, series="Naïve", PI=FALSE) +
  autolayer(snaive_fit_nc, series="Seasonal Naïve", PI=FALSE) + 
  autolayer(drift_fit_nc, series="Drift", PI=FALSE) +
  xlab("Year") + 
  ylab("Currency (€)") +
  ggtitle("Real Gross Domestic Product") +
  guides(colour=guide_legend(title="Forecast"))

rbind(
"Average" = accuracy(mean_fit_nc, test_nc)[2,c(2,3,5,6)],
"Naïve" = accuracy(naive_fit_nc, test_nc)[2,c(2,3,5,6)],
"Seasonal Naïve" = accuracy(snaive_fit_nc, test_nc)[2,c(2,3,5,6)],
"Drift" = accuracy(drift_fit_nc, test_nc)[2,c(2,3,5,6)])
##                     RMSE       MAE     MAPE      MASE
## Average        12157.667 10544.807 2.585153 1.5139682
## Naïve           7254.374  6161.831 1.528049 0.8846834
## Seasonal Naïve 10155.249  9703.442 2.403862 1.3931695
## Drift           6999.073  5140.116 1.289655 0.7379909

Il seasonal naïve method, pur essendo l’unico in grado di cogliere la stuttura stagionale della serie, presenta il secondo più elevato tasso di errore. Il più basso tasso di errore è fornito dal drift method, l’unico in grado di cogliere la componente di trend presente nei dati.

Natural cubic smoothing splines

Si può provare ad utilizzare un metodo di previsione più complesso e verificare se ciò porta ad un’accuratezza previsiva maggiore. Di seguito viene implementa la cubic smoothing spline, una forma di smoothing spline che utilizza funzioni cubiche a tratti per modellare la curva di regressione. Ogni osservazione rappresenta un nodo e, per evitare overfitting, vengono imposti vincoli sui coefficienti.

spline_fit_nc <- splinef(train_nc, 12)
autoplot(spline_fit_nc)

autoplot(gdp_nc) + 
  autolayer(mean_fit_nc, series="Mean", PI=FALSE) + 
  autolayer(naive_fit_nc, series="Naïve", PI=FALSE) +
  autolayer(snaive_fit_nc, series="Seasonal Naïve", PI=FALSE) + 
  autolayer(drift_fit_nc, series="Drift", PI=FALSE) +
  autolayer(spline_fit_nc, series="Spline", PI=FALSE) +
  xlab("Year") + 
  ylab("Currency (€)") +
  ggtitle("Real Gross Domestic Product") +
  guides(colour=guide_legend(title="Forecast"))

rbind(
"Average" = accuracy(mean_fit_nc, test_nc)[2,c(2,3,5,6)],
"Naïve" = accuracy(naive_fit_nc, test_nc)[2,c(2,3,5,6)],
"Seasonal Naïve" = accuracy(snaive_fit_nc, test_nc)[2,c(2,3,5,6)],
"Drift" = accuracy(drift_fit_nc, test_nc)[2,c(2,3,5,6)],
"Spline" = accuracy(spline_fit_nc, test_nc)[2,c(2,3,5,6)])
##                     RMSE       MAE     MAPE      MASE
## Average        12157.667 10544.807 2.585153 1.5139682
## Naïve           7254.374  6161.831 1.528049 0.8846834
## Seasonal Naïve 10155.249  9703.442 2.403862 1.3931695
## Drift           6999.073  5140.116 1.289655 0.7379909
## Spline          7830.295  6396.832 1.595106 0.9184236

Nonostante la complessità del modello ora utilizzato sia maggiore, le sue performance previsive non superano quelle di modelli semplici come il drift method o il naïve method.

Analisi dei residui

Per verificare se un modello ha catturato in modo adeguato le informazioni contenute nei dati, è necessario effettuare l’analisi dei residui, che dovrebbero avere le due seguenti proprietà: essere incorrelati, per mostrare che tutte le informazioni presenti nei dati sono state cattuarte dal modello, e avere media zero, per mostrare che le previsioni sono non distorte.
Inoltre, per la costruzione degli intervalli di confidenza, è utile che i residui abbiano varianza costante e siano normalmente distribuiti.
In sostanza, un buon modello è quello in cui i residui sono white noise.

Drift method

checkresiduals(drift_fit_nc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Random walk with drift
## Q* = 488.21, df = 8, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 8

Nonostante il metodo drift sia quello con una maggiore accuratezza previsiva, questo presenta residui altamente correlati e non è in grado di cogliere la stagionalità presente nei dati, visibile qui all’interno dei residui e del loro correlogramma. Come si poteva intuire, il test di Ljung-Box per l’incorrelazione dei coefficienti ha un p-value minore di \(\alpha=0.05\), quindi rifiutiamo l’ipotesi nulla di incorrelazione dei residui.
Inoltre, la distribuzione non è normale.

Seasonal Naïve Method

checkresiduals(snaive_fit_nc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 108.08, df = 8, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 8

Il metodo seasonal naïve, nonostante abbia una bassa accuratezza previsiva, presenta residui correlati solo per ritardi \(h=1, 2\) ed è in grado di cogliere la stagionalità presente nei dati. Inoltre, la sua distribuzione si può approssimare ad una normale. É presente asimmetria negativa a causa di due valori anomali. Infatti, anche dalla serie dei residui è evidente come il modello soffra della presenza degli outliers, legati alla crisi finanziaria del 2008. Anche in questo caso il test di Ljung-Box per l’incorrelazione dei coefficienti ha un p-value minore di \(\alpha=0.05\), quindi rifiutiamo l’ipotesi nulla di incorrelazione dei residui.

GDP in funzione di altre misure macroeconomiche

Per la costruzione di un modello per la previsione del GDP verranno di seguito prese in considerazione altre serie storiche macroeconomiche. Tuttavia, in questa sede non si procede con l’analisi dettagliata delle suddette serie (su cui ci sarebbe molto da dire), ma solo delle caratteristiche necessarie agli scopi dello studio.
Per ciascuna delle seguenti serie storiche, i dati estratti sono sullo stesso intervallo temporale considerato per la serie del GDP.

Real Exports of Goods and Services

Il Real Exports of Goods and Services si riferisce al valore totale di beni e servizi che l’Italia esporta verso altri Paesi, adeguato all’inflazione. Si tratta di un indicatore della competitività internazionale di un Paese.
Per l’Italia, le esportazioni sono una parte significativa dell’economia, poiché il Paese è noto per le sue esportazioni di prodotti di alta qualità.
Nota: Ci si riferirà a questa serie storica sempre con la sigla EXP.

Fonte: International Monetary Fund, Real Exports of Goods and Services for Italy [NXRNSAXDCITQ], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/NXRNSAXDCITQ, May 11, 2023.

exp <- ts(getSymbols("NXRNSAXDCITQ", src="FRED", auto.assign = FALSE), start=1996, frequency=4)
dyRangeSelector(dygraph(exp, 
                        main="Real Exports of Goods and Services",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Real Imports of Goods and Services for Italy

Il Real Imports of Goods and Services si riferisce al valore totale di beni e servizi che l’Italia importa da altri Paesi, adeguato all’inflazione. Si tratta di un indicatore economico della dipendenza dell’Italia dalle importazioni.
Come molte altre economie moderne, l’Italia importa una vasta gamma di beni e servizi.
Nota: Ci si riferirà a questa serie storica sempre con la sigla IMP.

Fonte: International Monetary Fund, Real Imports of Goods and Services for Italy [NMRNSAXDCITQ], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/NMRNSAXDCITQ, May 11, 2023.

imp <- ts(getSymbols("NMRNSAXDCITQ", src="FRED", auto.assign = FALSE), start=1996, frequency=4) 
dyRangeSelector(dygraph(imp, 
                        main="Real Imports of Goods and Services for Italy",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Early Estimate of Quarterly ULC Indicators: Total Labor Productivity

L’ULC, acronimo di Unit Labor Cost, è un indicatore macroeconomico che misura il costo del lavoro necessario per produrre un’unità di prodotto. L’ULC si calcola come il rapporto tra il costo del lavoro e la produttività del lavoro.
Questo indicatore può fornire informazioni importanti sulla capacità produttiva dell’economia italiana e sulla sua competitività internazionale.
Nota: Ci si riferirà a questa serie storica sempre con la sigla LAB.

Fonte: Organization for Economic Co-operation and Development, Early Estimate of Quarterly ULC Indicators: Total Labor Productivity for Italy [ULQELP01ITQ661N], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/ULQELP01ITQ661N, May 11, 2023.

lab <- ts(getSymbols("ULQELP01ITQ661N", src="FRED", auto.assign = FALSE) %>%
  window(start="1996-01-01"), start=1996, frequency=4) 
dyRangeSelector(dygraph(lab, 
                        main="Early Estimate of Quarterly ULC Indicators: Total Labor Productivity",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Real Final Consumption Expenditure for Italy

Con Real Final Consumption Expenditure ci si riferisce al valore totale della spesa finale effettuata dai residenti italiani per l’acquisto di beni e servizi, adeguato all’inflazione.
La spesa finale delle famiglie è un importante indicatore della domanda interna: una spesa finale elevata può indicare un aumento del reddito disponibile e una maggiore fiducia dei consumatori nell’economia. Tuttavia, una spesa finale eccessivamente elevata può anche indicare una maggiore inflazione o un aumento del debito delle famiglie.
Nota: Ci si riferirà a questa serie storica sempre con la sigla CONS.

Fonte: International Monetary Fund, Real Final Consumption Expenditure for Italy [NCRNSAXDCITQ], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/NCRNSAXDCITQ, May 11, 2023.

cons <- ts(getSymbols("NCRNSAXDCITQ", src="FRED", auto.assign = FALSE), start=1996, frequency=4) 
dyRangeSelector(dygraph(cons, 
                        main="Real Final Consumption Expenditure for Italy",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Real General Government Final Consumption Expenditure for Italy

Il Real General Government Final Consumption Expenditure si riferisce al valore totale della spesa finale effettuata dal governo per l’acquisto di beni e servizi (ad esempio le spese per la difesa, l’istruzione, la salute, la cultura, ecc.), adeguato all’inflazione.
La spesa pubblica può influire sull’economia nazionale attraverso il finanziamento di progetti infrastrutturali e di sviluppo, la creazione di posti di lavoro e l’aumento della domanda di beni e servizi. Tuttavia, una spesa pubblica eccessivamente elevata può anche portare a un aumento del debito pubblico e ad un aumento dei tassi di interesse.
Nota: Ci si riferirà a questa serie storica sempre con la sigla CONSGG.

Fonte: International Monetary Fund, Real General Government Final Consumption Expenditure for Italy [NCGGRNSAXDCITQ], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/NCGGRNSAXDCITQ, May 11, 2023.

consGG <- ts(getSymbols("NCGGRNSAXDCITQ", src="FRED", auto.assign = FALSE), start=1996, frequency=4)
dyRangeSelector(dygraph(consGG, 
                        main="Real General Government Final Consumption Expenditure for Italy",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Consumer Price Index: Total All Items for Italy

Il Consumer Price Index è un indicatore economico che misura la variazione dei prezzi di un paniere di beni e servizi consumati dalle famiglie. In pratica, l’indice dei prezzi al consumo fornisce una misura dell’inflazione in Italia.
L’aumento dei prezzi al consumo può avere un impatto significativo sulla capacità di acquisto delle famiglie e sull’andamento dell’economia.
Nota: Ci si riferirà a questa serie storica sempre con la sigla CPI.

Fonte: Organization for Economic Co-operation and Development, Consumer Price Index: Total All Items for Italy [CPALTT01ITQ657N], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/CPALTT01ITQ657N, May 10, 2023.

cpi <- ts(getSymbols("CPALTT01ITQ657N", src="FRED", auto.assign = FALSE) %>%
  window(start="1996-01-01", end="2022-10-01"), start=1996, frequency=4)
dyRangeSelector(dygraph(cpi, 
                        main="Consumer Price Index: Total All Items for Italy") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Real Gross Capital Formation for Italy

Con Real Gross Capital Formation ci si riferisce al valore totale degli investimenti effettuati dalle imprese, sia pubbliche che private, in beni di capitali al netto della deprezzamento e adeguato all’inflazione.
L’investimento in beni di capitale è una componente importante della spesa nazionale, poiché consente di aumentare la produttività e la capacità produttiva dell’economia italiana.
Nota: Ci si riferirà a questa serie storica sempre con la sigla CAP.

Fonte: International Monetary Fund, Real Gross Capital Formation for Italy [NIRNSAXDCITQ], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/NIRNSAXDCITQ, May 11, 2023.

cap <- ts(getSymbols("NIRNSAXDCITQ", src="FRED", auto.assign = FALSE), start=1996, frequency=4)
dyRangeSelector(dygraph(cap, 
                        main="Real Gross Capital Formation for Italy",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Costruzione del dataset

data <- cbind(gdp, exp, imp, lab, cons, consGG, cpi, cap)
colnames(data) <- c("gdp", "exp", "imp", "lab", "cons", "consGG", "cpi", "cap")

Per evitare i problemi legati alla multicollinearità, bisogna assicurarsi che le variabili non siano altamente correlate tra loro.

ggpairs(as.data.frame(data))

Risulta opportuno escludere le variabili imp e consGG, rispettivamente altamente correlate alle variabili exp e cons.

data <- cbind(gdp, exp, lab, cons, cpi, cap)
colnames(data) <- c("gdp", "exp", "lab", "cons", "cpi", "cap")
ggpairs(as.data.frame(data))

Dagli scatterplot risulta evidente l’esistenza di una relazione lineare tra la variabile GDP e le covariate prese in considerazione, suggerendo di poter costruire un modello di regressione lineare.

autoplot(data, facets = TRUE)

Modello di regressione lineare multiplo

Il modello più semplice che possiamo considerare è il modello di regressione lineare, che assume la presenza di un legame lineare tra la variabile gdp e le covariate sopra presentate. Inoltre, si aggiungono le variabili trend e season, legate alla struttura della serie storica dipendente.

\[ gdp = \beta_0 + \beta_1 exp + \beta_2 lab + \beta_3 cons + \beta_4 cpi + \beta_5 cap + \beta_6 trend + \beta_7 season + \epsilon \]

mod <- tslm(gdp ~ exp + lab + cons + cpi + cap + trend + season, data=data)
autoplot(gdp, series="Data") +
  autolayer(fitted(mod), series="Fitted") +
  ylab("GDP (€)")

Il modello sembra adattarsi molto bene ai dati, ma è necessario analizzare le sue statistiche di sintesi.

summary(mod)
## 
## Call:
## tslm(formula = gdp ~ exp + lab + cons + cpi + cap + trend + season, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4612.2 -1128.3   -94.5  1330.7  3438.6 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.443e+04  1.187e+04   2.058 0.042259 *  
## exp          4.583e-01  4.093e-02  11.196  < 2e-16 ***
## lab          4.127e+02  1.559e+02   2.647 0.009452 ** 
## cons         7.741e-01  2.479e-02  31.222  < 2e-16 ***
## cpi         -4.871e+02  3.709e+02  -1.313 0.192145    
## cap          5.126e-01  3.528e-02  14.532  < 2e-16 ***
## trend       -1.196e+02  3.050e+01  -3.922 0.000163 ***
## season2      5.416e+03  5.898e+02   9.182 7.21e-15 ***
## season3      4.230e+03  6.477e+02   6.531 2.91e-09 ***
## season4      5.847e+03  7.474e+02   7.822 6.07e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1891 on 98 degrees of freedom
## Multiple R-squared:  0.991,  Adjusted R-squared:  0.9902 
## F-statistic:  1202 on 9 and 98 DF,  p-value: < 2.2e-16

Un dato particolarmente rivelante è che tutte le variabili, ad eccezione del cpi (indice dei prezzi al consumo) sono statisticamente significative. Ciò è confermato anche dal fatto che il p-value della F-statistic ci porta a rifiutare l’ipotesi nulla \(H_0\) di nullità di tutti i predittori e ad accettare l’ipotesi alternativa \(H_1\) di almeno un predittore non nullo.

Il modello così costruito presenta un elevatissimo indice \(R^2\) corretto, una misura della bontà di adattamento. L’indice \(R^2\) indica la percentuale di varianza nella variabile dipendente spiegata dalle variabili indipendenti, ma non tiene conto del numero di variabili indipendenti utilizzate nel modello, il che può portare a una sovrastima dell’efficacia del modello quando vengono utilizzate molte variabili indipendenti. L’indice \(R^2\) corretto aumenta solo se l’aggiunta di una variabile indipendente migliora effettivamente la capacità del modello di spiegare la variazione nella variabile dipendente.
L’indice \(R^2\) corretto, così come quello standard, ha un valore compreso tra 0 e 1: quanto più il valore è prossimo a 1, tanto migliore è l’adattamento del modello ai dati. In questo caso l’indice corretto \(R^2= 0.99\) mostra un ottimo adattamento del modello ai dati.

Tuttavia, è indispensabile effettuare un’analisi dei residui. Dalle statistiche sopra riportate possiamo dedurre che la distribuzione dei residui è leggermente asimmetrica positiva e presenta uno standard error, misura della varianza dei residui, molto elevato. Si deduce che il modello di regressione non è in grado di spiegare tutta la variabilità presente nella variabile dipendente.

Analisi dei residui

Per i motivi sopra menzionati, è importante effettuare un’analisi dei residui per valutare il modello di regressione.

checkresiduals(mod)

## 
##  Breusch-Godfrey test for serial correlation of order up to 13
## 
## data:  Residuals from Linear regression model
## LM test = 58.143, df = 13, p-value = 1.126e-07

I residui mostrano eteroschedasticità, appaiono correlati e il p-value del test di Breusch-Godfrey è molto minore di 0.05, quindi non c’è evidenza empirica per poter affermare che i residui sono in linea con un processo white noise. Tuttavia, la correlazione sembra essere forte solo per i ritardi \(h=1,2\) e la distribuzione si può approssimare ad una normale.

Modello di regressione lineare a tratti

Si potrebbe ipotizzare che questo problema sia legato ai cambi di pendenza nella serie storica del GDP, dunque si procede con la definizione di un modello lineare a tratti.

dyRangeSelector(dygraph(gdp, 
                        main="Real Gross Domestic Product",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )

Si procede con l’utilizzo di due nodi, uno in corrispondenza del quarto trimestre del 2007 e l’altro del primo trimestre del 2014.

t <- time(gdp)
tau.2007 <- c(2007,4)
tau.2014 <- c(2014,1)
tau1 <- ts(pmax(0, t - tau.2007), start = 1996)
tau2 <- ts(pmax(0, t - tau.2014), start = 1996)
mod_pw <- tslm(gdp ~ t + tau1 + tau2 + exp + lab + cons + cpi + cap + season)
summary(mod_pw)
## 
## Call:
## tslm(formula = gdp ~ t + tau1 + tau2 + exp + lab + cons + cpi + 
##     cap + season)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4244.5  -970.6   -86.9  1152.4  3952.4 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.475e+06  3.712e+05   3.973 0.000137 ***
## t           -7.121e+02  1.826e+02  -3.901 0.000178 ***
## tau1         1.836e+02  2.393e+02   0.768 0.444650    
## tau2        -3.777e+02  1.476e+02  -2.558 0.012090 *  
## exp          5.526e-01  4.418e-02  12.509  < 2e-16 ***
## lab         -1.806e+01  1.806e+02  -0.100 0.920590    
## cons         7.951e-01  2.621e-02  30.332  < 2e-16 ***
## cpi         -2.896e+02  3.477e+02  -0.833 0.406889    
## cap          5.238e-01  3.361e-02  15.582  < 2e-16 ***
## season2      3.963e+05  2.535e+05   1.563 0.121256    
## season3      4.090e+03  6.027e+02   6.786 9.48e-10 ***
## season4      3.977e+05  2.536e+05   1.568 0.120129    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1756 on 96 degrees of freedom
## Multiple R-squared:  0.9924, Adjusted R-squared:  0.9915 
## F-statistic:  1142 on 11 and 96 DF,  p-value: < 2.2e-16
checkresiduals(mod_pw)

## 
##  Breusch-Godfrey test for serial correlation of order up to 15
## 
## data:  Residuals from Linear regression model
## LM test = 43.828, df = 15, p-value = 0.000117

Il modello presenta ancora i residui correlati, ma il p-value del test di Breusch-Godfrey è aumentato, pur restando sotto il livello \(\alpha=0.05\). Si potrebbe, pertanto, pensare di aggiungere altri due noti, consapevoli del rischio di overfitting.

dyRangeSelector(dygraph(gdp_clean, 
                        main="Real Gross Domestic Product",
                        ylab="Currency (€)") %>%
                  dyEvent("2020-04-01", "Covid", labelLoc = "bottom") %>%
                  dyEvent("2009-01-01", "Crisi", labelLoc = "bottom")
                )
t <- time(gdp)
tau.2007 <- c(2007,4)
tau.2014 <- c(2014,1)
tau.2018 <- c(2018,4)
tau.2020 <- c(2020,3)
tau1 <- ts(pmax(0, t - tau.2007), start = 1996)
tau2 <- ts(pmax(0, t - tau.2014), start = 1996)
tau3 <- ts(pmax(0, t - tau.2018), start = 1996)
tau4 <- ts(pmax(0, t - tau.2020), start = 1996)
mod_pw2 <- tslm(gdp ~ t + tau1 + tau2 + tau3 + tau2 + exp + lab + cons + cpi + cap + season)
summary(mod_pw2)
## 
## Call:
## tslm(formula = gdp ~ t + tau1 + tau2 + tau3 + tau2 + exp + lab + 
##     cons + cpi + cap + season)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4252.4  -966.2   -77.6  1126.5  4025.0 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.498e+06  3.727e+05   4.020 0.000116 ***
## t           -7.238e+02  1.833e+02  -3.949 0.000151 ***
## tau1         3.008e+02  2.754e+02   1.092 0.277495    
## tau2        -6.924e+02  3.936e+02  -1.759 0.081773 .  
## tau3         2.251e+02  2.609e+02   0.863 0.390465    
## exp          5.533e-01  4.425e-02  12.505  < 2e-16 ***
## lab         -9.249e+00  1.812e+02  -0.051 0.959391    
## cons         7.882e-01  2.746e-02  28.707  < 2e-16 ***
## cpi         -3.299e+02  3.513e+02  -0.939 0.350055    
## cap          5.359e-01  3.646e-02  14.696  < 2e-16 ***
## season2      3.421e+05  2.615e+05   1.308 0.193975    
## season3      4.200e+03  6.169e+02   6.809 8.84e-10 ***
## season4      3.435e+05  2.616e+05   1.313 0.192409    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1759 on 95 degrees of freedom
## Multiple R-squared:  0.9925, Adjusted R-squared:  0.9915 
## F-statistic:  1044 on 12 and 95 DF,  p-value: < 2.2e-16
checkresiduals(mod_pw2)

## 
##  Breusch-Godfrey test for serial correlation of order up to 16
## 
## data:  Residuals from Linear regression model
## LM test = 43.329, df = 16, p-value = 0.0002492

Il modello presenta ancora i residui correlati, ma il p-value del test di Breusch-Godfrey è aumentato, pur restando sotto il livello \(\alpha=0.05\).

Bisogna fare un’importante precisazione: la correlazione nei residui non rende le previsioni distorte, quindi queste non sono sbagliate, rende solo gli intervalli di previsione più ampi e quindi la stima più incerta.

rbind("lineare multiplo" = CV(mod),
      "lineare a tratti v.1" = CV(mod_pw),
      "lineare a tratti v.2" = CV(mod_pw2))
##                           CV      AIC     AICc      BIC     AdjR2
## lineare multiplo     4203785 1641.201 1643.951 1670.705 0.9901983
## lineare a tratti v.1 3757857 1627.032 1630.904 1661.899 0.9915435
## lineare a tratti v.2 3799479 1628.189 1632.705 1665.739 0.9915209

Siccome in termini di accuratezza previsiva tutti e tre i modelli presentati offrono ottime performance, il modello prescelto è il più semplice tra i due modelli lineari a tratti.

Plot dei residui

Per assicurarsi che il modello lineare sia adatto, si può verificare che lo scatter dei residui non mostri alcun pattern.

df <- as.data.frame(data)
df[,"Residuals"]  <- as.numeric(residuals(mod_pw))
p1 <- ggplot(df, aes(x=exp, y=Residuals)) +
  geom_point()
p2 <- ggplot(df, aes(x=lab, y=Residuals)) +
  geom_point()
p3 <- ggplot(df, aes(x=cons, y=Residuals)) +
  geom_point()
p4 <- ggplot(df, aes(x=cpi, y=Residuals)) +
  geom_point()
gridExtra::grid.arrange(p1, p2, p3, p4, nrow=2)

Al netto di alcuni valori anomali, i residui sembrano essere randomly scattered.

Previsioni

Su test set

É possibile utilizzare il modello di regressione lineare multiplo costruito mod per fare previsioni e valutare la bontà previsiva del modello. Per i motivi di cui sopra, il modello viene addestrato sulle osservazioni dal 1996 al 2016 e testato sulle osservazioni dal 2017 al 2019 della serie storica originale gdp.

nc_train <- ts(data[1:84,], start=1996, frequency = 4)
nc_test <- as.data.frame(ts(data[85:96,], start=2017, frequency = 4))
mod_nc <- tslm(gdp ~ exp + lab + cons + cpi + cap + trend + season, data=nc_train)
fcast_mod <- forecast(mod_nc, nc_test)
autoplot(gdp_nc) + 
  autolayer(fcast_mod, series="Forecast", PI=FALSE) +
  xlab("Year") + 
  ylab("Currency (€)") +
  ggtitle("Real Gross Domestic Product")

accuracy(fcast_mod, test_nc)[2,c(2,3,5,6)]
##         RMSE          MAE         MAPE         MASE 
## 2736.1123039 2360.1509437    0.5810829    0.3388581

Sul futuro

Possiamo, inoltre, utilizzare il modello così costruito per fare previsione sui valori futuri di GDP, previsioni che non possiamo testare.

fcast_mod_new <- forecast(mod$fitted.values, h=12)
autoplot(fcast_mod_new) +
  ggtitle("Forecasting") +
  xlab("Year") +
  ylab("Currency (€)") 

Grazie!

Appendice

Di seguito vengono riportati alcuni dati relativi ai decessi causati dal Covid-19.
I dati sono solo stati solo visulizzati a scopo esplorativo al fine di avere una visione migliore su come il Covid abbia potuto influenzare il GDP.
Non è stata fatta modellazione con questi dati in quanto aggregando la serie da settimanale a trimestrale, come la serie GDP oggetto di studio, le osservazioni erano così poche da non permettere di fare inferenza.

covid <- read.csv("/Users/mariapiatedesco/Downloads/HEALTH_MORTALITY_11052023211912210.csv") %>%
  select(c(WEEK, YEAR, Value)) %>%
  arrange(YEAR, WEEK)
covid_ts <- ts(covid[,3], start=c(2020,5), frequency = 52)
autoplot(covid_ts)

ggsubseriesplot(covid_ts)

covid_new <- covid %>%
  add_row(WEEK=1, YEAR=2020, Value=0) %>%
  add_row(WEEK=2, YEAR=2020, Value=0) %>%
  add_row(WEEK=3, YEAR=2020, Value=0) %>%
  add_row(WEEK=4, YEAR=2020, Value=0) %>%
  filter(YEAR<2023) %>%
  arrange(YEAR, WEEK) %>%
  select(Value)
covid_new_ts <- ts(covid_new, start=c(2020,1), frequency = 52)
covid_new_q <- aggregate(covid_new_ts, nfrequency = 4)
autoplot(covid_new_q)

ggsubseriesplot(covid_new_q)

gdp_covid <- ts(gdp[89:100], start=c(2020,1), frequency = 4)
covid_data <- cbind(covid_new_q, gdp_covid)
autoplot(covid_data, facets=TRUE)

cor.test(gdp_covid, covid_new_q)
## 
##  Pearson's product-moment correlation
## 
## data:  gdp_covid and covid_new_q
## t = 0.44283, df = 10, p-value = 0.6673
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4728543  0.6600501
## sample estimates:
##       cor 
## 0.1386816